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CALCULATION  OF  A SUPERSONIC  FLOW 
AROUND  A COOLABLE  SPHERE  USING  THE 
COMPLETE  NAVI ER-STOKES  EQUATIONS 

B.  M.  Pavlov 


This  work  [1]  presents  a statement  and  a numerical  solution  for 
a problem  of  a stationary  flow  around  blunt  solids  by  a super- 
sonic flow  of  a low-density  gas.  The  method  is  based  on  a finite- 
difference  approximation  of  the  complete  nonstationary  Navier-Stokes 
equations  for  a compressible  gas  and  on  finding  a stationary  solu- 
tion as  a result  of  the  determination.  Works  [1]— [3]  present  the 
initial  calculation  results  and  examine  certain  methodological 
problems  which  characterize  the  numerical  algorithm.  Work  [4] 
contains  the  calculation  results  for  a flow  past  a heat-insulated 
sphere.  This  work  discusses  a pattern  of  flow  around  a coolable 
sphere  under  the  conditions  of  small  Reynolds  numbers. 

Designations:  r,  0 - spherical  polar  system  of  coordinates 

(Fig.  1);  r,  - sphere  radius;  W - flow  velocity  vector;  u,  v - its 
w 

components  with  respect  to  r and  0;  p,  h,  T,  p - density,  enthalphy , 
temperature,  and  pressure;  A,  y - coefficients  of  heat  conductivi- 
ty and  shear  viscocity;  y - ratio  of  specific  heats;  q - heat 
flow  to  the  wall;  t - friction  tension;  F - force  of  resistance; 

CD  - coefficient  of  total  resistance;  C„  - friction  coefficient; 


1 - mean  free  path  of  molecules;  Fie,  W»,  Kn,  Pr  - numbers  of 
Reynolds,  Mach,  Knudsen,  and  Prandl ; d,  I - thickness  and  intensi- 
ty of  a shock  wave;  pc  - stagnation  pressure  behind  a normal  shock; 
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rA  - sonic  line;  subscripts:  00  - in  an  incident  flow,  w - on  the 
sphere  surface,  s - behind  a normal  inviscid  shock,  0 - flow  stag- 
nated adiabatically  and  isentropically , nc  - in  a boundary  layer, 

CM  - in  a free-molecular  flow;  superscripts:  0 - in  the  leading 
critical  point,  * - along  the  boundary  ray  0*. 

1.  An  examination  is  made  of  a viscous  gas  flow,  which  is  uniform 
and  supersonic  at  infinity,  around  a sphere.  The  problem  concerning 
the  calculation  of  the  flow  field  is  solved  under  the  following 
assumptions:  1)  the  gas  flowing  past  the  sphere  is  thermically  and 
calorically  ideal,  2)  the  coefficients  \ and  y are  power  functions 
of  enthalpy  (~A“,  0,5  *£  u>  «£  1),  3)  the  flow  has  a laminar  nature  at 
small  Reynolds  numbers,  and  4)  under  these  conditions  of  flow  the 
Navier-Stokes  equations  for  the  compressible  gas  are  valid  in  the 
entire  area  outside  the  sphere. 

Conditions  1)  and  2)  lead  to  constant  values  of  the  numbers  y 

and  Pr.  Dimensionless  variables  are  used:  linear  dimensions  per- 

2 

tain  to  rw,  velocities  - to  Woo,  density  -^to  p<» , enthalpy  - to  , 
pressure  and  tension  of  friction  - to  P00W00 , heat  flow  - to  p^Woo , 
and  force  of  resistance  - to  l/2jx/£p,.lFi. 

Work  [1]  shows  the  initial  system  of  equations  in  a dimension- 
less form  relative  to  the  unknowns  p,  u,  v,  h,  p,  and  y in  the 
system  of  coordinates  r,  9.  Since  with  a laminar  nature  the  flow 
has  a symmetry  relative  to  the  9=0  axis,  it  is  sufficient  to  seek 
a solution  for  the  problem  in  the  region  r^l , 040-$tt.  According  to 
the  systematic  investigations,  which  are  described  briefly  in  [1] 
and  [3],  numerical  integration  can  be  accomplished  in  the  end  region 
(Fig.  1)  l^r4R2j  O^0£0*,  where  R2=R2(6)  - given  eauation  for  the 
surface  of  rotation  (in  this  case,  o^  an  ellipsoid  whose  left  focal 
point  is  congruent  with  the  center  of  the  sphere),  0*  - position 
of  the  boundary  ray  (0<9*,$tt). 

Let  us  formulate  the  boundary  conditions  of  the  problem.  The 

conditions  of  the  uniformity  of  the  flow  at  an  infinite  distance 

from  the  sphere  are  removed  to  the  surface  R2(0),  O-S04P*.  On  the 

sphere  the  enthalpy  (temperature)  of  the  surface  coincides  with 

the  enthalpy  of  gas  and  is  equal  to  the  value  h (there  is  no  jump 

w 

in  temperature)  and  both  velocity  components  are  equal  to  zero 
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(there  is  no  slip).  On  the  axis  0=0,  B=n  we  use  the  conditions  of 
symmetry  of  the  flow.  If  9*<n,  it  is  necessary  to  assign  conditions 
for  the  functions  u,  v,  h along  the  ray  0=0*.  By  means  of  series 
of  calculations  with  a varied  position  of  the  angle  6* ,O<0*<tt  it 
was  established  that,  when  assigning  the  approximate  conditions 
for  these  functions  in  the  form  of  o2 f/O02=O  (i.  e.,  a linear  change 
in  f in  the  vicinity  of  the  ray  0*)  there  is  an  insignificant  reverse 
effect  of  the  flow.  In  this  case  the  density  p(r,  0*)  can  be  found 
from  an  analogous  condition  or  from  a continuity  equation. 

The  finite-difference  approximation  of  the  initial  eouations  and 
boundary  conditions,  and  also  the  peculiarities  of  the  numerical 
algorithm,  are  described  in  [1]— [33* 

2.  The  study  of  the  super-  and  hypersonic  flow  around  blunt  solids 
with  the  cooling  of  their  surface  (h  <h  ) is  of  practical  interest 

W O 

in  connection  with  the  problem  of  heat  protection  for  a flying 

vehicle.  Below  we’ll  discuss  the  characteristic  peculiarities  of 

the  viscous  gas  flow  near  the  sphere  placed  in  a supersonic  flow 

of  a diatomic  ideal  gas  (y  = l,  4 Pr=0.72,  co=3/4,  h =0.3).  The  ef- 

w 

feet  of  the  Moo  and  Reoo  numbers  has  been  studied  on  the  flow  pattern 
under  the  conditions  of  flow  and  heat  exchange,  which  are  indicated 
in  the  table.  The  same  table  shows  the  Knudsen  Kn  and  Reynolds  Res 
(calculated  using  the  parameters  behind  the  normal  shock)  numbers 
which  characterize  the  degree  of  raref cation  of  the  incident  flow: 


Rem  _ Rem 

Hi  ~ 


[ Mv  + » 1» 

[ 2y-  (V-1)/*L  J ' 


With  an  increase  in  rarefaction  there  is  an  increase  in  the  Kn 
number  and  a decrease  in  Res<  In  most  calculations  0*=99°. 

Let’s  briefly  describe  the  behavior  of  the  gas-dynamic  parameters 
at  the  leading  critical  point  in  the  area  where  the  shock  wave  is 
thickened  and  in  a compressed  layer  behind  it,  and  also  along  the 
surface  of  the  sphere.  Since  the  calculations  were  accomplished 
for  the  flow  conditions  which  occupy  the  intermediate  position 
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between  the  flows  of  the  solid  and  free-molecular  media.  It  Is  of 
interest  to  compare  the  results  of  the  calculation  with  the  corres- 
ponding theoretical  results  and  experiments  in  rarefied  gases. 

1.  Leading  critical  point.  The  solid  line  in  Fig.  2,  a repre- 
sents the  behavior  of  the  ratio  p°/p',  which  characterizes  the 

w o 

deviation  of  the  calculated  pressure  from  the  pressure  in  an  in- 
viscid  flow,  as  a function  of  Re,  vr.  . With  a decrease  in  the 

parameter  n = Re,Vp,  (introduced  for  a convenience  of  the  compari- 

son with  the  experimental  data  [5],  shaded  area  in  this  figure)  the 
values  of  p°/p ' begin  to  increase  approximately  from  n=100,  which 
is  in  good  agreement  with  the  experiments  and  the  theory  of  free- 
molecular  flows  [6]  (with  Re  =0,  y = 1.4  and  totally  diffused  deflec- 
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tion  of  the  molecules  from  a hard,  very  cold  wall  the  ratio  pCM/p^ 
attains  a value  equalling  to  1.087).  An  increase  in  pressure  is 
the  result  of  merging  of  the  thickened  boundary  layer  and  shock 
wave . 

The  solid  line  in  Fig.  2,  b shows  the  behavior  of  the  ratio 
q°/q°  as  a function  of  the  Re  number.  This  ratio  characterizes 
the  difference  in  heat  transfer  in  a flow  of  rarefied  gas  and  in 
flows  with  large  Reynolds  numbers.  The  value  q°  was  calculated 
by  the  Fay  and  Riddell  formula  [7],  and 

«.<»>-- i^rO*  •£-).• 

According  to  the  theory  of  free-molecular  flow  q°^=l/2  and  the 

ra  0 ^CM/c^nc,  which  is  proportional  to  V Re.,  approaches  zero  with 
an  increase  in  rarefaction.  The  same  figure  shows  the  following: 
the  shaded  area  - data  of  a series  of  experiments  (borrowed  from 
[5],  the  dot-dash  line  - q°^/q°c , the  dashed  line  - the  curve 

&/ <ftic=  I +0,52 /V Re„  obtained  in  [8]  by  processing  the  measurements 

(the  conditions  of  this  experiment  are  close  to  the  conditions  of 
these  calculations:  Moo  *2 , 4,  6;  20£Res£800  , Tw/T  «0.  33) . In  most 
of  the  experiments  carried  out  to  determine  heat  transfer  the 

value  q°/q°  has  the  following  nature  of  change:  with  a decrease 

w nc 

in  Re.  10»>  this  value  increases,  here  <><?nc  , attains 
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a maximum  at  certain  Re,  102)  , and  then  decreases 

(q£<Vnc  when  /te,<:iO}.  The  calculated  values  of  <£/ ^ detect  a 

similar  behavior;  although  they  begin  to  decrease  when  He  *10,  they 

s 

still  remain  higher  than  those  determined  in  the  experiments  (with 

the  same  Re  number).  Apparently,  at  low  Reynolds  numbers.  Re  <102 , 
s s 

one  should  take  into  account  in  the  calculation  the  conditions  of 

slip  and  temperature  jump  on  the  wall,  which,  in  particular,  lead 

to  a decrease  in  the  values  of  the  heat  flow  [9]. 

With  a decrease  in  the  Re„  number  the  density  p°  increases  in 

o ^ ^ 

proportion  to  p since  the  h is  constant  (with  thermal  insulation 
w w 

a drop  was  observed  in  the  density  at  the  leading  critical  point 
[4]  and  [10]). 

2.  Shock  wave  structure.  The  thickness  of  a shock  wave  in  a 
viscous  gas  is  usually  determined  by  the  formula  proposed  by  Prandl : 

d h-tm 

(df!dr)  mu 


where  f - any  of  the  gas-dynamic  functions,  df/dr  - maximum  slope 
of  the  profile  f to  the  axis  0=0.  In  most  of  the  experiments  car- 
ried out  to  study  the  shock-wave  structure  (for  example,  [11]  and 
[12])  the  density  field  was  measured  and,  therefore,  thickness  d 
was  calculated  using  the  profile  p in  the  area  of  the  shock  wave. 
However,  under  the  cooling  conditions  of  the  surface  the  calculated 
density  profiles  p(r,  0)  in  the  shock  wave  have  a complex  behavior 
(see  Figs.  4 and  5)  and  it  is  difficult  to  use  them  for  the  calcula- 


tion of  d.  For  this  it  is  more  convenient  to  use  the  profiles 
u(r,  0)  and  p(r,  0),  which  have  approximately  the  same  nature  as 
with  thermal  insulation.  The  calculation  of  the  thickness  o^  the 
shock  wave  using  the  profiles  u(r,  0),  when  woo=6  and  Reoo=300,  200, 
90,  45,  and  20,  yield  the  following  values:  d=0.055;  0.084;  0.169; 
0.286;  and  0.474.  Thus,  with  an  increase  in  rarefaction  the  shock 
wave  thickens  quickly.  The  comparison  of  the  thickness  d in  cases 
of  thermal  insulation  and  cooling  (with  all  other  conditions  being 
equal)  indicates  a somewhat  lesser  thickness  in  the  case  of  cooling. 

The  thickness  of  the  transition  zone  in  the  shock  wave,  when  there 
is  a change  in  the  angle  0 , can  be  determined  conditionally  by  the 


value  D(0  )=r-^(0)— r2(0  ) > where  r^  and  r^  are  points  at  the  fixed 
ray  0:  in  r^  the  value  p— p^xO.Ol  p^ , in  r2  the  dpjdr  attains  a 
maximum  value  in  the  shock  wave.  The  position  of  the  washed  out 
shock  wave  relative  to  the  sphere  can  then  be  determined  by  the 
middle  line  r^  of  the  indicated  zone.  The  dot-dash  curve  on  Fig. 

6 shows  the  position  of  the  line  r^  when  Re0O=45  and  1^=6.  This 
line  travels  in  approximately  the  same  way  as  a shock  wave  during 
an  inviscid  flow  but  it  is  located  further  from  the  sphere.  A 
similar  behavior  of  the  shock  wave  zone  was  observed  in  an  ex- 
periment (see  [ 5 ] ) • With  a change  in  the  numbers  and  Re^  the 
behavior  of  the  value  A=r^— 1,  which  can  be  roughly  considered  as 
the  "deoarture"  of  the  shock  wave,  has  the  following  feature: 

a)  with  a decrease  in  the  number  Reoo  (Moo,0  are  constant)  A increases 

b)  with  an  increase  in  0 (Re^,  are  constant)  A increases;  c) 
with  a change  in  the  number  M*,  (PjR^  are  constant)  this  function 
is  not  monotonic  (see  Fig.  4):  for  2^M„<4  the  A diminishes,  for 
Moo>4  it  increases  (the  same  behavior  was  observed  in  the  calcula- 
tions of  a thermoinsulated  sphere  [4]).  The  shock  wave  intensity 

f 

was  determined  by  the  value  of  the  ratio  l(0)=(p  — p«>)/(P  — p ), 

max  o 

where  p is  the  value  of  maximum  pressure  in  the  area  of  the 
max 

shock  wave  at  the  fixed  ray  0.  It  was  found  that  with  an  increase 
in  0 the  function  1(0),  for  all  the  flow  conditions  examined,  de- 
creases (for  example,  when  Moo=6 ; /(0)  =0,917;  / ^-2- j =0,238  for 

Rea>  = 200  and  /(0)  =0,882,  / fy)  =0,193  forReco=45). 

Thus,  with  an  increase  in  the  angle  0 and  a decrease  in  the  Re^, 
number,  there  is  a continuous  (according  to  r)  washing  out  and 
weakening  of  the  shock  wave  (just  as  in  the  case  of  thermal  insula- 
tion [ l 3— [ 33). 

3.  Structure  of  a comnressed  layer.  Figure  3 (a,  b,  c)  shows 
a typical  behavior  of  the  flow  parameters  in  a compressed  layer 
(Meo=8,  Reoo=90),  where,  for  comparison,  the  dashed  lines  show  the 
solution  when  0=0  in  the  case  of  thermal  insulation  of  the  sohere. 
From  an  analysis  of  the  calculated  data  we  draw  the  following  con- 
clusions. The  enthalpy  h(r,  0)  has  a maximum  on  each  rav  0 due  to 
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the  cooling  of  the  sphere's  surface.  The  density  p(r,  9)  has  a 

local  maximum  immediately  behind  the  shock  wave  as  it  moves  away 

from  the  region  of  the  leading  critical  point,  and  maximum  pressure 

p(r,  0)  is  attained  between  the  shock  wave  and  wall;  with  the  low 

values  of  0 the  maximum  pressure  p(r,  0)  is  on  the  wall.  With  a 

relatively  large  Re  numbers  (for  example,  when  M00=2  and  Reoo=500) 

s 

there  is  an  area  close  to  the  sphere's  surface  in  which  the  pres- 
sure along  the  normal  to  the  wall  changes  insignificantly,  which 
indicated  the  existance  of  a region  of  the  type  of  a boundary  layer. 

The  position  of  the  sonic  line  r#(0)  is  depicted  in  Fig.  6 for 
the  various  Re^  (1^=6,  solid  lines)  numbers  and  in  Fig.  8 for  the 
various  (Reoo=90,  solid  lines)  numbers.  It  is  evident  that  with 
a change  in  0 and  (Rem  is  constant)  the  behavior  of  these  lines 
is  nonmonotonic.  For  a constant  Pico  number  sonic  lines  move  further 
away  from  the  sphere's  surface  with  a decrease  in  the  Re<„  number. 

Let  us  examine  in  more  detail  the  behavior  of  the  solution  along 
the  stagnation  line  0=0.  Figure  4 characterizes  the  behavior  of 
u and  p with  a change  in  the  (Reoo=90)  number.  Using  profiles 
u(r,  0)  we  can  clearly  see  the  nonmonotonic  nature  of  "departure" 
of  the  shock  wave  noted  above.  The  density  increases  along  the 
0 = 0 axis  when  r -*•  1;  here  there  are  regions  of  sharp  increase  of 
p in  the  shock  wave  and  close  to  the  body  (the  latter  is  the  result 
of  cooling  of  the  surface).  Both  these  regions  have  a tendency  to 
merge  with  an  increase  in  the  M^,  number  and  the  fixed  Reoo  number 
(i.  e.,  with  a decrease  in  the  Reg  number)  (Figs.  4 and  5). 

Figure  5 shows  the  distribution  of  the  density  p(r,  0)  and 

enthalpy  h(r,  0)  when  Me =6  and  a change  in  the  Re^  number.  It  is 

evident  that  with  a decrease  of  the  Re«>  number  (and,  consequently. 

Re  number)  the  maximum  value  of  h(r,  0),  which  is  achieved  behind 

the  shock  wave,  also  decreases  and  h <h  . If  we  take  the  point 

max  s 

r,  where  the  maximum  h(r,  0)  is  achieved,  to  be  the  trailing  front 
of  the  washed-out  shock  wave,  it  turns  out  that,  at  this  point,  the 
density  p(r,  0)<p  and  p(r,  0)<p  , i.  e.,  an  incomplete  comnres- 
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sion  in  the  shock  wave,  leads  to  disruption  of  the  Renkin-Hugoniot 
conditions  and  to  its  gradual  disintegration  with  an  increase  in 
rarefaction.  With  an  increase  in  rarefaction  its  gradual  thicken- 

j 


ing  is  evident  in  the  area  of  the  shock  wave  according  to  the 
profiles  of  h(r,  0).  Approximately  the  same  behavior  of  the  density 
and  tempera' ure  was  observed  along  the  line  of  stagnation  in  the 
experiment  of  work  [13],  in  which  nitrogen  flowed  past  a cylinder 
with  a spherical  blunt  end  (M00=9,  Kn=0.045,  To=298°K,  and  l^/T^O.3). 

The  field  temperatures  determined  in  [3]  attest  to  the  presence  of 

a noticeable  temperature  jump  on  the  model  surface  (T°/T  =0.78  when 

w o 

T /T  «0.3  at  the  leading  critical  point), 
w o 

4.  Distributions  along  the  surface  of  the  snhere.  Figure  7 

shows  the  graphs  of  the  deduced  pressure  Pv=pv(B)lpv( 0)  and  heat 

flow  q (0-)  when  M =6  and  with  different  values  of  the  Rea,.  In  the 
w 

same  figure  the  dot-dash  line  depicts  the  function  n^^=cos0/2,  the 

dotted  line  - <?nc(0),  obtained  from  the  calculation  of  the  boundary 

layer  on  the  sphere  [14]  (^*=10,  fle~106,  y=  1.4,  Tw/T0=0,6),  the 

dashed  line  - Pu,(B)  when  Mco=6,  Reoc  = <x>  ([7],  chapter  8).  With 

a decrease  in  the  Rea,  number  the  values  of  p gradually  increase 

on  the  front  section  of  the  sphere  and  they  are  higher  than  the 

corresponding  values  in  an  inviscid  flow.  With  an  increase  in  the 

angle  0 the  function  q (0)  decreases  gradually  and  with  a decrease 

w 

in  the  Re^,  the  values  of  q (0)  increase;  in  this  case,  <7nc<<7«><<7cM. 

w 

The  nature  of  change  in  the  calculated  heat  flow  is  in  agreement 
with  the  experiments  [8]  and  [15]- 

Figure  8 shows  the  friction  coefficient  graphs  Cf  = tb / p^wL  = 

Re^'  (n  dv/dr)",  for  Reco=90  and  the  various  Moo  numbers  (solid  lines). 

In  the  same  figure  the  dot-dash  line  depicts  the  function 
[C/(0)]cM  = (sin20)/2.  the  dashed  line  - [C/(0)]nc  taken  from 

work  [14].  All  calculation  curves  of  C^(0)  are  located  between  the 
curves  for  the  free-molecular  and  boundary-layer  conditions  of 
flow.  The  values  of  increase  with  a decrease  in  the  Reoo  number. 
With  a change  in  the  angle  0 the  function  C^.(0)  has  a maximum  for 
the  various  Re,*,  numbers  when  0-40—60  degrees,  i.  e.,  about  in  the 
middle  of  the  area  of  a sharp  approch  of  the  sonic  line  to  the 
sphere.  In  the  calculations  of  the  boundary  layer  on  blunted  solids 
[16]  it  was  found  that  the  maximum  values  of  the  function  C^(9) 
were  obtained  in  the  vicinity  of  the  sonic  point  on  the  surface 
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of  the  streamlined  body.  The  behavior  of  the  functions  p and  C„ 

*w  f 

is  analogous  to  the  case  of  thermal  insulation  when  there  is  a 
decrease  in  the  Re^  number  [4]. 

Figure  9 shows  the  construction  of  the  ratio  0^/(0^) C!VI  as  a 

function  of  the  number  2Re  , the  shaded  areas  indicate  the  data  of 

s 

the  experiments  carried  out  to  determine  the  resistance  of  spheres 

e* 

[51.  The  quantity  CD  was  calculated  by  the  formula  CD  = 4 j (pa cos 6 + 

0 

Cf  sin  0) sin 6 <f6,  and  (Cd)qM  by  the  formula  used  to  calculate  the 
diffused  deflection  with  the  accommodation  coefficient  eaualling 
unity  [6].  It  is  obvious  that  in  the  examined  conditions  of  flow 
Cd<(Cd)cm  ’ while  the  ratio  being  examined  increases  with  a de- 
crease in  the  Re  number. 

s 

2 

3.  The  numerical  results  indicate  that  wner:  Re  .£10  the  shock 

s 

wave  and  the  boundary  layer  join  in  the  area  of  the  leading  critical 
point  of  the  sphere,  which  leads  to  a series  of  new  phenomena  which 
do  not  take  place  in  flows  with  large  Reynolds  numbers.  With 
reduced  rarefaction  the  calculated  aerodynamic  characteristics 
have  a tendency  to  approach  the  corresponding  values  in  a free- 
molecular  flow  and,  in  the  range  10-^Res^102J  there  is  a satisfac- 
tory agreement  between  them  and  the  measurements  derived  from  the 
experiments.  When  Re  s-$  10  there  is  a certain  anomaly  in  the  behavior 
of  some  characteristics  of  the  flow,  which  can  be  explained  by  the 
fact  that  the  effects  of  slip  and  temperature  jump  were  disregarded 
in  these  calculations . Apparently,  when  taking  these  effects  into 
consideration  within  the  framework  of  the  Navier-Stokes  model,  there 
should  be  a better  agreement  with  the  experiments  also  at  lower 

Reynolds  numbers  (Re  <10 ).  The  same  conclusions  are  d- a.m  from  the 

s 

analysis  of  flows  with  thermal  insulation  of  the  sphere's  surface 
[4].  Thus,  one  can  assert  that  the  numerical  solutions  of  the 
problem  of  flow  in  this  setting  afford  the  description  of  the 
phenomena  in  a compressed  layer,  which  are  close  to  those  observed 
in  experiments. 

The  author  expresses  his  appreciation  to  Z.  M.  Yemel’yanov  and 
L.  V.  Kuznetsov  who  took  part  in  calculations  and  processing  of  the 
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results . 
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